A WEAK GALERKIN FINITE ELEMENT METHOD WITH 
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Abstract. The novel idea of weak Galcrkin (WG) finite element methods is on the use of weak 

functions and their weak derivatives defined as distributions. Weak functions and weak derivatives 

can be approximated by polynomials with various degrees. Different combination of polynomial 

spaces leads to different weak Galerkin finite element methods, which makes WG methods highly 

ff^ fiexible and efficient in practical computation. This paper explores the possibility of optimal com- 

^_^ bination of polynomial spaces that minimize the number of unknowns in the numerical scheme, yet 

^> without compromising the accuracy of the numerical approximation. For illustrative purpose, the 

^vj authors use second order elliptic problems to demonstrate the basic idea of polynomial reduction. A 

new weak Galerkin finite element method is proposed and analyzed. This new finite element scheme 

features piecewise polynomials of degree fc > 1 on each element plus piecewise polynomials of degree 

fc — 1 > on the edge or face of each element. Error estimates of optimal order are established 

■^r for the corresponding WG approximations in both a discrete H^ norm and the standard L^ norm. 

In addition, the paper presents a great deal of numerical experiments to demonstrate the power of 

\l the WG method in dealing with finite element partitions consisting of arbitrary polygons in two di- 

Cm mensional spaces or polyhedra in three dimensional spaces. The numerical examples include various 

finite element partitions such as triangular mesh, quadrilateral mesh, honey comb mesh in 2d and 

' J ' mesh with deformed cubes in 3d. The numerical results show a great promise of the robustness, 

■^T reliability, flexibility and accuracy of the WG method. 
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1. Introduction. This paper is concerned with weak Galerkin (WG) finite cl- 

1— H ement methods by exploring optimal use of polynomial approximating spaces. In 

^ general, weak Galerkin refers to finite element techniques for partial differential equa- 

" tions in which differential operators (e.g., gradient, divergence, curl, Laplacian) are 

I approximated by weak forms as distributions. The main idea of weak Galerkin finite 

w~) element methods is the use of weak functions and their corresponding discrete weak 

• derivatives in algorithm design. For the second order elliptic equation, weak func- 

/-^ tions have the form of w = {vo,vi,} with v — vq inside of each element and v = Vb 

^^ on the boundary of the element. Both vq and Vf, can be approximated by polynomi- 

1-H als in PiiT) and Ps{e) respectively, where T stands for an element and e the edge 

►^ or face of T, £ and s are non-negative integers with possibly different values. Weak 

• ^ derivatives are defined for weak functions in the sense of distributions. For com- 

r% puting purpose, one needs to approximate the weak derivatives by polynomials. For 

2 example, for the weak gradient operator, one may approximate it in the polynomial 

space [Prn{T)]'^. Various combination of (P^(T),Ps(e), [Pm(T)]'') leads to different 

class of weak Galerkin methods tailored for specific partial differential equations. The 
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goal of this paper is to explore optimal combination of the polynomial spaces Pe(T) 
and Psie) that minimizes the number of unknowns without compromising the rate of 
convergence for the corresponding WG method. 

For simplicity, we demonstrate the idea of optimality for polynomials by using 
the second order elliptic problem that seeks an unknown function u satisfying 

(1.1) - V • (aVw) = /, in f7, 

(1.2) u — g, on dil, 

where f2 is a polytopal domain in M'' (polygonal or polyhedral domain for d = 2,3), 
Vu denotes the gradient of the function u, and a is a symmetric d x d matrix- valued 
function in fi. We shall assume that there exists a positive number A > such that 

(1.3) Ca^ > AC*C, Ve e R'^. 

Here ^ is understood as a column vector and ^* is the transpose of ^. 

A weak Galerkin method has been introduced and analyzed in [S] for second 
order elliptic equations based on a discrete weak gradient arising from local RT [^ or 
BDM jT] elements. More specifically, in the case of BDM element of order fc > 1, the 
gradient space is taken as [Pm{T)]'^ = [Pk{T)]'^ and the weak functions are defined by 
using {Pi{T),Ps{e)) = (Pfe-i(r), Pfc(e)). For the i?T element of fc > 0, the gradient 
space is the usual RT element for the vector component while the weak functions 
are given by (Pf(T),P^(e)) = {Pk{T) , Pk{e)) . Due to the use of the RT and BDM 
elements, the WG finite element formulation of [8 is limited to classical finite element 
partitions of triangles {d = 2) or tetrahedra {d = 3). In addition, the corresponding 
WG scheme exhibits a close connection with the standard mixed finite element method 



for (1.1|-(1.2) 



The main goal of this paper is to investigate the possibility of optimal combination 
of polynomial spaces that minimize the number of unknowns in the numerical scheme 
without compromising the order of convergence. The new WG scheme will use the 
configuration of (Pfe(T), Pk-i{ e) , P k -iiT Y) , and the corresponding WG solution con- 



verges to the exact solution of (1.1)-(1.2) with rate oiO{h'') in H^ and 0(/i'^+^) in L"^ 
norm, provided that the exact solution of the original problem is sufficiently smooth. 
It should be pointed out that the unknown vq associated with the interior of each 
element can be eliminated in terms of the unknown Vb defined on the element bound- 
ary in practical implementation. This means that, for problems in K^, only edges of 
the finite element partition shall contribute unknowns (fc unknowns from each edge) 
to the global stiffness matrix problem. The new WG scheme is, therefore, a natural 
extension of the classical Crouzix-Raviart Pi non-conforming triangular element to 
arbitrary order and arbitrary polygonal partitions. 

It have been proved rigorously in [3] that Pk type of polynomials can be used 
in weak Galerkin finite element procedures on any polygonal/polyhedral elements. It 
contrasts to the use of polynomials Pk for triangular elements and tensor products Qk 
for quadrilateral elements in classic finite element methods. In practice, allowing ar- 
bitrary shape in finite element partition provides a great flexibility in both numerical 
approximation and mesh generation, especially in regions where the domain geometry 
is complex. Such a flexibility is also very much appreciated in adaptive mesh refine- 
ment methods. Another objective of this paper is to study the reliability, flexibility 
and accuracy of the weak Galerkin method through extensive numerical tests. The 
first and second order weak Galerkin elements are tested on partitions with differ- 
ent shape of polygons and polyhedra. Our numerical results show optimal order of 



convergence for fc = 1,2 on triangular, quadrilateral, honey comb meshes in 2d and 
deformed cube in 3d. 

One close relative of the WG finite element method of this paper is the hybridiz- 
able discontinuous Galerkin (HDG) method 4 . But these two methods are funda- 
mentally different in concept and formulation. The HDG method is formulated by 
using the standard mixed method approach for the usual system of first order equa- 
tions, while the key to WG is the use of discrete weak differential operators. For the 



second order elliptic problem (1.1 )-( 1.2 ), these two methods share the same feature of 
approximating first order derivatives or fluxes through a formula that was commonly 
employed in the mixed finite element method. For high order PDEs, such as the 
biharmonic equation [6], the WG method is greatly different from the HDG. It should 
be emphasized that the concept of weak derivatives makes WG a widely applicable 
numerical technique for a large variety of partial differential equations which we shall 
report in forthcoming papers. 

The paper is organized as follows. In Section [2l we shall review the definition of 
the weak gradient operator and its discrete analogues. In Section [3] we shall describe 
a new WG scheme. Section ID will be devoted to a discussion of mass conservation for 
the WG scheme. In Section [5] we shall present some technical estimates for the usual 
L^ projection operators. Section [6] is used to derive an optimal order error estimate 
for the WG approximation in both H^ and L^ norms. Finally in Section [7J we shall 
present some numerical results that confirm the theory developed in earlier sections. 

2. Weak Gradient and Discrete Weak Gradient. Let K be any polytopal 
domain with boundary dK. A weak Junction on the region K refers to a function 
V = {vq^vi,} such that Vq € Li^iK) and vi, £ H^[dK). The first component Vq can be 
understood as the value of v in K, and the second component Vb represents ■;; on the 
boundary of K. Note that Vb may not necessarily be related to the trace of vq on dK 
should a trace be well-defined. Denote by W{K) the space of weak functions on K; 
i.e., 

(2.1) WiK) ^{v = {vo^Vb} : v^ e L\K), Vb G H^OK)}. 

Define {v,w)d = jj^vwdx and {v,w) — J vwds. 

The weak gradient operator, as was introduced in [H], is defined as follows for the 
completion of the paper. 

Definition 2.1. The dual of L^{K) can be identified with itself by using the 
standard L^ inner product as the action of linear functionals. With a similar inter- 
pretation, for any v G W{K), the weak gradient ofv is defined as a linear functional 
V^v in the dual space of H{div,K) whose action on each q S H{div,K) is given by 

(2.2) (V^w,g)K = -(wo,V ■q)K + {vb,q-n)gK, 

where n is the outward normal direction to dK , (uq, V • q)K = /i^ i'o(V • q)dK is the 
action of Vq onV ■ q, and {vb,q ■ vl)qk is the action of q ■ n on Vb £ H^i^dK). 

The Sobolev space H^{K) can be embedded into the space W{K) by an inclusion 
map iw '■ H^{K) — > W{K) defined as follows 

iw{'P) = {<f>\K,cl)\aK}, cj,eH\K). 

With the help of the inclusion map iw, the Sobolev space H^{K) can be viewed as a 
subspace of W{K) by identifying each (f> e H^{K) with iw{4>)- Analogously, a weak 



function v = {vo,Vb} £ W{K) is said to be in H^{K) if it can be identified with 
a function (j) G H^{K) through the above inclusion map. It is not hard to see that 
the weak gradient is identical with the strong gradient (i.e., Vu,w = Vf) for smooth 
functions v G H^{K). 

Denote by Pr{K) the set of polynomials on K with degree no more than r. We 
can define a discrete weak gradient operator by approximating V^, in a polynomial 
subspace of the dual of H{div, K). 

Definition 2.2. The discrete weak gradient operator, denoted by '^w.r.K, is de- 
fined as the unique polynomial {^ w,r.Kv) G [Pr{K)Y satisfying the following equation 

(2.3) {V^^r,KV,q)K^-{vo,V ■q)K + {vb,q-Xl)aK, 'i q C, [Pr{K)Y ■ 

By applying the usual integration by part to the first term on the right hand side 



of (2.3 1, we can rewrite the equation (2.3 1 as follows 



(2.4) {Vw,r,KV, q)K = (Vuo, q)K + {vb ~Vo,q- n)oK, Vg G [Pr{K)f- 

3. Weak Galerkin Finite Element Schemes. Let 7/i be a partition of the 
domain fl consisting of polygons in two dimension or polyhedra in three dimension 
satisfying a set of conditions specified in [S] . Denote by £/i the set of all edges or flat 
faces in Th, and let f ° = £h\dfl be the set of all interior edges or flat faces. For every 
element T E Th, we denote by hr its diameter and mesh size h = max^gT;, ft-T for Th- 

For a given integer fc > 1, let V/i be the weak Galerkin finite element space 
associated with Th defined as follows 

(3.1) Vh = {v = {vo,Vb} : vo\t e Pk{T), Vb\e e Pk-iie), eedT,Te %} 
and 

(3.2) Vl^ ^{v: veVh, vb = on dn}. 

We would like to emphasize that any function v ^ Vh has a single value Vb on each 
edge ee£h- 

For each element T E Th, denote by Qq the L^ projection from L'^{T) to Pk{T) 
and by Qb the L^ projection from L^(e) to Pk-i{e). Denote by Qh the L^ projection 
onto the local discrete gradient space [Pfe_i(T)]''. Let V — H^{il). We define a 
projection operator Qh ■ V -^ Vh so that on each element T ETh 

(3.3) QhV ^ {QoVo,QbVb}, {v(i,Vb} = iw{v) ^W{T). 

Denote by '^ w.k-i the discrete weak gradient operator on the finite element space 



Vh computed by using (2.3 1 on each element T; i.e. 



(V^,,fe_lu)|T = ^iu,k-l,T{v\T), Vw e Vh- 

For simplicity of notation, from now on we shall drop the subscript /c — 1 in the 
notation 'Ww,k-i for the discrete weak gradient. 
Now we introduce two forms on Vh as follows: 

a(v,w) == ^ (aVu,w, V^,w)t, 
TeTh 

s{v,w) ^ p ^ h^^{QbVo - Vb, QbWo - Wb)dT, 

TGTh 



where p can be any positive number. In practical computation, one might set p — 1- 
Denote by as(-, •) a stabihzation of a(-, •) given by 

as{v, w) = a{v, w) + s{v, w). 



Weak Galerkin Algorithm 1. A numerical approximation for (1.1) and (1.2) 
can be obtained by seeking ut — {uojUfe} G Vh satisfying both Uf, ~ Qi,g on dfl and 
the following equation: 



(3.4) 



,{uh,v) = {f,vo), y v = {vo,vb} ev°. 



Note that the system (3.4 1 is symmetric and positive definite for any parameter 
value of p > 0. 



Next, we justify the well-postedness of the scheme (3.4). For any v £ Vh, let 

(3.5) nielli := ^/as{v,v). 

It is not hard to see that ||| • ||| defines a semi-norm in the finite element space Vh. We 
claim that this semi-norm becomes to be a full norm in the finite element space V^. 
It suffices to check the positivity property for ||| • |||. To this end, assume that v G V^ 
and |||u||| = 0. It follows that 

(aV^w,Vu,i;) +P X] ^T^iQbVo - Vb,QbVo ~ Vb)aT = 0, 
TeTh 

which implies that V^w = on each element T and QbVQ = Vb on dT. It follows from 
V„w = and K^ that for any q e [Pk-iiT)]'^ 



= {VwV,q)T 
= (Vwo, <?)t - {vo -Vb,q- n)oT 
= (Vwo, q)T - {QbVQ -Vb,q- ujqt 
= (Vwo,<7)t- 

Letting q ~ Vuo in the equation above yields Vwq = on T G 7h- Thus, uo ~ const 
on every T G T/i- This, together with the fact that Qb^o — Vb on dT and f f, = on 
Oil, implies that vq = Vb = 0. 



tion 



Lemma 3.1. The weak Galerkin finite element scheme (3.4) has a unique solu- 
Proof If u)j and uf^ are two solutions of (3.4), then eh — Uh — u^ would 



satisfy the following equation 

a,ieh,v)^0, yveV,^. 
Note that eh G V^. Then by letting u = e^ in the above equation we arrive at 



\\eh\\\ = as(e/i,e/i) = 0. 



It follows that e/i = 0, or equivalently, Uf^ 
lemma. D 



(1) - .,(2) 



This completes the proof of the 



4. Mass Conservation. The second order elliptic equation ( |1.1[ ) can be rewrit- 
ten in a conservative form as follows: 



V • g = /, q 



-aVu. 



Let T be any control volume. Integrating the first equation over T yields the following 
integral form of mass conservation: 



(4.1) 



/ qnds= I fdT. 

JdT Jt 



We claim that the numerical approximation from the weak Galerkin finite element 



method (3.4) for (1.1 1 retains the mass conservation property (4.1) with an appropri 



ately defined numerical flux q^. To this end, for any given T S Th, we chose in (3.4 1 



a test function v — {vq, Vb = 0} so that z^o = 1 on T and Vq = elsewhere. It follows 
from (lOI that 



(4.2) 



aS/^Uh ■ \7wvdT + phj, / {Qtuo — ui,)ds 



fdT. 



Let Qh be the local L^ projection onto the gradient space [Pk-i{T)]'^. Using the 
definition (2.31 for V^w one arrives at 



Jt Jt 



(4.3) 



}h{aVwUh)dT 



/ QhiaWwUfi) ■ nds. 

JdT 



Substituting ( |4.3[ ) into ( |4.2[ ) yields 

(4.4) / {-Qh (aW^Uh) + ph;^\QbUn - Uh)n} ■ nds ^ f fdT, 

JdT Jt 

which indicates that the weak Galerkin method conserves mass with a numerical flux 
given by 



qh 



h (aV^.Wh) + phj}{Q},UQ - Ub)n. 



Next, we verify that the normal component of the numerical flux, namely qh • n, 
is continuous across the edge of each element T. To this end, let e be an interior 
edge/face shared by two elements Ti and T2. Choose a test function v = {vo,Vb} so 
that "^o = and vj, = everywhere except on e. It follows from (3.4) that 



(4.5) 



/ aVwUh-^wvdT-phrj} I {QbUo - Ub)\TiVbds 

JT1UT2 JdTine 



= 



T2 / {QbUQ - Ub)\T2Vbds 

dT2ne 



Using the definition of weak gradient ( |2.3[ ) we obtain 

/ aVyjUh ■ VwvdT = / Qh(aV„,u/i) • V ^vdT 

JTiUTs JT1UT2 

= / iQhiaVwUh)\Ti ■ni+Qh{aV^Uh)\T2 ■n2)vbds, 

J e 

wliere n^ is ttie outward normai direction of Ti on the edge e. Note that rii + n2 
Substituting the above equation into (|4.5| yields 



/ (-Q/i(aV^u/i)|Ti + phj}{QbUo - Wfc)|Tini) • niv^ds 

J e 

(-Q;i(aV^M;J|T2 + pKj^HQbUQ - Ub)|T2n2) • n2Wfcds, 



which shows the continuity of the numerical flux qh in the normal direction. 

5. Some Technical Estimates. This section shall present some technical re- 
sults useful for the forthcoming error analysis. The first one is a trace inequality 
established in [5] for functions on general shape regular partitions. More precisely, let 
T be an element with e as an edge. For any function </? e H^{T), the following trace 
inequality holds true (see [S] for details): 

(5.1) Ml<C{h^'M'T + hT\\\7ipfT)- 

Another useful result is a commutativity property for some projection operators. 

Lemma 5.1. Let Qh and Q/j be the L^ projection operators defined in previous 
sections. Then, on each element T (z Th, we have the following commutative property 

(5.2) V^(Q,,0)=Q,,(V0), y^eH^iT). 



Proof. Using (2.3), the integration by parts and the definitions of Qh and Qh, we 



have that for any r € [Pk^i{T)]'^ 

(V^(Q,,0), r)T = -{Qo(f>, W-t)t + {Qb<l>, T ■ n)oT 
= -(0, \/ ■t)t + (0, r • n)QT 
= (V(^,t)t 
= (Q„(V(/.),t)t, 



which implies the desired identity (5.2). D 



The following lemma provides some estimates for the projection operators Qh 
and Qh. Observe that the underlying mesh Th is assumed to be sufficiently general 
to allow polygons or polyhedra. A proof of the lemma can be found in [9] . It should 
be pointed out that the proof of the lemma requires some non-trivial technical tools 
in analysis, which have also been established in [5]. 

Lemma 5.2. LetTk be a finite element partition offl that is shape regular. Then, 
for any (j) G H''~^^[n), we have 

(5.3) E ll<^ - Qo<^IIt + E ^tIIVI^/- - Qo^T < Ch'^'^+^'^UWl^,, 

Ten Ten 

(5.4) J2 ii«(^'^ - 'Qhi^miT < ch^'ml+i- 

Ten 



Here and in what follows of this paper, C denotes a generic constant independent of 
the meshsize h and the functions in the estimates. 

In the finite element space Vh, we introduce a discrete H^ semi- norm as follows: 

(5.5) \\v\\i,h^{^ {\\'^vofT + h-j,^\\QbVt)~Vb\\dT)\ ■ 

\TeTh J 



The following lemma indicates that || • ||i_/,, is equivalent to the trip-bar norm (3.5 1. 

Lemma 5.3. There exist two positive constants Ci and C2 such that for any 
V = {vQ,Vb\ € Vh, we have 

(5.6) Ci\\v\\i^H<\M<C2\\v\\i.h. 



Proof. For any v = {ug, Vb\ G V/i, it follows from the definition of weak gradient 



(2.4) and Qb that 



(5.7) (V^t;,g)T = (Vfo,g)T + (w6-QbWo,g-n)aj,, \/q e[Pk-i{T)Y. 



By letting q = Vu,u in (5.7 1 we arrive at 



(Vu,u, V^,v)t = (V-wo, Vi„u)t + {vb - QbVQ, '^wV ■ n)gj,. 



From the trace inequality (5.1) and the inverse inequality we have 



(V^w, Vu,u)t < ||Vwo||t||V^w||t + WQbVa - VbWarW'^wvWdT 

< ||VwoI1t||V^,w||t + C/i^^/^llQbWo - Vb\\aT\\ywv\\T 



Thus, 

||V„i;|1t < C iWVvoWl + hr'WQbVo - VblUr) ^ , 
which verifies the upper bound of |||w|||. As to the lower bound, we chose q = Vuq in 



(5.7) to obtain 

(W^v, Vwo)t = (Vwq, Vuo)t + {vb - QbVQ, Vuo • n)^^,. 
Thus, from the trace an inverse inequality we have 

llVwolll < IIV^wIItIIVvoIIt + Ch-^/^\\QbVo - VfcHaTllV^ollT. 
This leads to 

||Vt;o||T < C (||V^«||| + Ch^'WQbVo - Wb|||r)' , 

which verifies the lower bound for |||w|||. Collectively, they complete the proof of the 
lemma. D 

Lemma 5.4. Assume that Th is shape regular. Then for any w G _ff''+^(J7) and 
V = {vo,Vb} G Vh, we have 

(5.8) \siQhW,v)\ < Ch''\\w\\k+i\\\vl 

(5.9) \£Uv)\<Ch''\\w\\k+i\\\v\\\, 



where £w{v) = Y^reTn (a(^^ ~ QftVw) • n, uq - Wf,)aT- 



Proof. Using the definition of Qb, (|5.1|), and (|5.3|), we obtain 

\s{QhW,v)\ 



y^ hj}{Qb{Qow) - Qbw, Q 



bVa ~ Vb)dT 



Ten 



^ hrj}{QaW-w, QbVo-Vb)dT 
Ten 



<ciY. (hr^WQow - w\\j. + llV(Qo^ - u.)|||) 
XTeTh > 

I X! hj.^\\QbVo-Vb\\9T I 
VTeTh / 

< Ch''\\w\\k+i\\\v\\\. 



As to (5.9), it follows from the Cauchy-Schwarz inequality, the trace inequality 



(5.1) and the estimate (5.4) that 



(5.10) \e^{v)\ = 



^ (a(Vu; - QftVw) • n, wq - Vb)aT 



Ten 



Ten 

< c [ ^ hrMVw ~ Q,yw)\\lT] I E ^T^ho - ybWdT 

\Ten / \Ten J 

yren ) 



Using the trace inequality (5.1) and the approximation property of the L^ projection 
operator we obtain 

||wo - I'bllaT < Iko - QfaWollaT + IIQbWo - I'bllaT 

1 /2 

<Chrj{ \\\7vo\\T + \\QbVQ-Vb\\aT- 



Substituting the above inequality into (5.101 yields 



(5.11) K{v)\ < Ch'^Mk+i i J2 {W^^oWT + hT^WQbVo-VbWlr} 

\Ten / 



which, along with the estimate (5.6), verifies the desired estimate (5.9). D 



6. Error Analysis. The goal of this section is to establish some error estimates 



for the weak Galerkin finite element solution Uh arising from (3.4). The error will 



be measured in two natural norms: the triple-bar norm as defined in (3.5) and the 
standard L^ norm. The triple bar norm is essentially a discrete H^ norm for the 
underlying weak function. 



10 



For simplicity of analysis, we assume that the coefficient tensor a in ( |1.1[ ) is a 
piecewise constant matrix with respect to the finite element partition 7/i. The result 
can be extended to variable tensors without any difficulty, provided that the tensor a 
is piecewise sufficiently smooth. 

6.1. Error equation. Let Uh = {uo,Ub} € Vh be the weak Galerkin finite 



element solution arising from the numerical scheme (3.4 1. Assume that the exact 
solution of ( 1.1 1-( 1.2 ) is given by u. The L^ projection of u in the finite element space 



Vh is given by 



Let 



QhU = {Qou,Qbu}. 



eh = {eo, ej = {Qqu - uq, Q^u - uj 

be the error between the WG finite element solution and the L^ projection of the 
exact solution. 

Lemma 6.1. Let e^ be the error of the weak Galerkin finite element solution 



arising from {3.4-)- Then, for any v G Vf^ we have 

(6.1) as{eh,v) =£^{v) + s{QhU,v), 

where £u{v) = J^TeT^ (a(Vu - Q/jVu) • n, wo - Vb)9T- 

Proof. Testing (|l.l[) by using uq of u = {vq, Vb} G V^ we arrive at 



(6.2) 



^ (aV-u, Wvq)t - ^ (aVw • n, uq - Vb)dT = (/, vq), 
Ten Ten 



where we have used the fact that J^TeT {^'^'^ ' ^i'^b)dT = 0. To deal with the term 
J^TeT ('^^'"i Vwo)t in (6.2 1, we need the following equation. For any (j) e H^{T) and 
V E Vh, it follows from (U^l, the definition of the discrete weak gradient (2.3|, and 



the integration by parts that 

{aW^Qh(f>,^wv)T = (aQ/i(V0),Vt„i;)T 

= -(fo, V • (aQftV0))T + {vb, (aQ,,V0) • n)aT 
= (V?;o,aQftV0)T - («o - Vb, (aQ/iV0) • n)oT 
(6.3) = (aV(/), Vwo)t - ((aQ/iV(/)) • n, vq - Wb)^. 

By letting (f> = u in (6.3), we have from combining (6.3) and (6.2) that 



^ {aVyjQhU, Vyjv)T = (/, Vq) + ^ (a(Vu - Q^Vm) • n, uo - Ub)aT 
TeTh TeTh 

= (/,«o)+4(i^). 

Adding s{QhU, v) to both sides of the above equation gives 

(6.4) OsiQhU, v) = (/, Wo) + tu{v) + s{QhU, v). 



Subtracting (3.4 1 from (6.4) yields the following error equation, 

as{eh,v) ^iuiv) + s{QhU,v), V-w G V^°. 



This completes the proof of the lemma. D 
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6.2. Error estimates. The error equation (6.1) can be used to derive the fol- 



lowing error estimate for the WG finite element solution. 

The orem 6.2 . Let Uh G Vh be th e weak Galerkin finite element solution of the 
problem (1.1)-\1.2) arising from (3.4)- Assume the exact solution u G H^'^^[U,). 



Then, there exists a constant C such that 



(6.5) 



\\uh-Qhu\\ < Ch''\\u\\k+i. 



Proof. By letting v = e^ in (|6.1[), we have 



(6.6) 



||e/i||| ^ luieh) + siQhU, e/j). 



It then follows from (5.8) and (5.9) that 



||eh||| < C/i'' II u|lfc+i III e;i III, 



which implies (6.5). This completes the proof. D 



Next, we will measure the difference between u and Uh in the discrete H^ semi- 
norm II • Ijiji as defined in (5.5). Note that (5.5) can be easily extended to functions 



in iJ^(il) + Vfi through the inclusion map iw- 

Coroll ary 6 .3. Let Uh S Vh be th e weak Galerkin finite element solution of 
the problem \1.1)-(1.2) arising from (3.4)- Assume the exact solution u G _ff'^+^(r2). 
Then, there exists a constant C such that 



(6.7) 



\u-Uh\\i,h<Ch''\\u\\ 



k+l- 



Proof It follows from (5.6) and (6.5) that 



WQhU - Uh\\i,h < C\\\QhU - u/ill < Ch''\\u\\k+i- 
Using the triangle inequality, ( |5.3[ ) and the equation above, we have 

\\u-Uh\\i,h < \\u- Qhu\\i.h + WQhU- UhWiji < Ch'^\\u\\k+i. 
This completes the proof. D 

In the rest of the section, we shall derive an optimal order error estimate for the 
weak Galerkin finite element scheme (3.4) in the usual L^ norm by using a duality 



argument as was commonly employed in the standard Galerkin finite element methods 
[3JI2]. To this end, we consider a dual problem that seeks $ G -f^o(^) satisfying 



(6.8) 



V • (aV$) ^ eo in fl. 



Assume that the above dual problem has the usual _ff^ -regularity. This means that 
there exists a constant C such that 



(6.9) 



I<i'll2 < qjeoll. 



The orem 6.4 . Let Uh G Vh be t he weak Galerkin finite element solution of the 
problem (1.1)-(1.2) arising from {3. 4)- Assume the exact solution u G H^^^{n). In 
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addition, assume that the dual problem (6.8) has the usual H^ -regularity . Then, there 
exists a constant C such that 



(6.10) 



\u~uo\\<Ch^+^\\u\\k+i. 



Proof. By testing (6.8) with eg we obtain 

lleof = -(V-(aV<i>),eo) 



(6.11) 



= ^ (aV$, Veo)T " ^ (aV$ • n, eo - e,,)^^,, 
TeTh TeTh 



where we have used the fact that Cf, = on dVi. Setting = <J> and v = Ch in (6.31 
yields 

(6.12) (aV^Q,,$, V„eft)T - (aV$, Veo)T - ((aQ/.V$) • n, eo - efc)^^,. 



Substituting (6.121 into (6.11) gives 



(aV^.e,,, V^Q,,$) + ^ (a(Q,,V$ - V$) • n, eo - eb)^^, 
(6.13) =(aV^,e;„ V^Q,,$)+^$(eft). 



|eo| 



It follows from the error equation ( |6.1[ ) that 

(6.14) {aV^eu, V^Qh^) = iuiQu"^) + s{QhU, Qh^) - s{eh, Qh^)- 



By combining (6.13) with (6.14) we arrive at 

(6.15) ||eo||2 = iuiQh'^) + s{Qhu, Qh'^) - s{eh, Q^*) + ^*(e,0. 



Let us bound the terms on the right hand side of (6.151 one by one. Using the 
triangle inequality, we obtain 



KiQum 



< 



(6.16) 



Y^ {a{Vu - Qh^u) ■ n, Qo* - Qb^) 

""en 

J2 (a(Vu - QhVu) • n, Qo* - <^)ot 
Y^ {a{Vu - QhVu) • n, $ - Qb^)9T 



dT 



Ten 



We first use the definition of Qb and the fact that $ = on dfl to obtain 
(6.17) J2 (a(Vu-QhVii)-n,$-Qb$)aT== ^ (aVw • n, $ - Qb$)gy = 0. 



Ten 



Ten 



From the trace inequality (5.1) and the estimate (5.3) we have 

1/2 

- 



5]||Qo$-$|||t <ch^m 

\Ten / 
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and 



[ ^ ||a(Vu - Q„Vu)||2^ ) < Ch''---\\u\\k+i. 
\TeTh / 

Thus, it follows from the Cauchy-Schwarz inequality and the above two estimates that 



J2 (a(Vw - Qh^u) ■ n, Qo$ - $)aT 



TeTh 



XTeTh J yreTh J 



(6.18) < Ch^+'Wullk+iim 



Combining (6.16) with (6.171 and (6.18) yields 

(6.19) K{Qh'^)\ < Ch''+'\\u\\k+i\m\ 



Analogously, it follows from the definition of Qb, the trace inequality (5.1), and the 



estimate (5.3 1 that 



\s{QhU, (9/,.$)| < p^ hrj.^ \{Qb{Qou) - Qbu, QbiQo'^) - Qb^)aT\ 

Ten 

<pY1 hr'WQbiQou - u)\\9T\\QbiQo'^ - '^)\\dT 

Ten 
<pY1 K^WQou-uWarWQo'^-'^WaT 



TGTh 



1/2 



VTGTh / \TeTh J 



1/2 



(6.20) 



<c/i*=+i||u|U+i||$|| 



The estimates (5.8) with fc = 1 and the error estimate (6.5) imply 
(6.21) \s{eh. g„$)| < C/i||$|i2|||e„||| < Ch!'^^\\u\\k^^\\<^h. 



Similarly, it follows from (5.9) and (6.5) that 

(6.22) K$(e,)|<C;^'=+l|l^i|U+l|l<I>|l, 



Now substituting ( |6.19[ )-( |6.22| into ( |6.15| yields 



which, combined with the regularity assumption (6.9 1 and the triangle inequality 



gives the desired optimal order error estimate (6.10). D 
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7. Numerical Examples. In this section, we examine the WG method by test- 
ing its convergence and flexibihty for solving second order elliptic problems. In the 
test of convergence, the first (fc = 1) and second (fc = 2) order of weak Galerkin 
elements are used in the construction of the finite element space Vh- In the test of 
flexibility of the WG method, elliptic problems are solved on finite element partitions 
with various configurations, including triangular mesh, deformed rectangular mesh, 
and honeycomb mesh in two dimensions and deformed cubic mesh in three dimensions. 
Our numerical results confirm the theory developed in previous sections; namely, op- 
timal rate of convergence in H^ and L^ norms. In addition, it shows a great flexibility 
of the WG method with respect to the shape of finite clement partitions. 

Let Uh — {uo,Ub} and u be the solution to the weak Galerkin equation and the 
original equation, respectively. The error is defined by e^ = u^ — QhU = {eo,eb}, 
where eg = Uq ~ Qo"^ a-nd ei, = ui, — QbU- Here QhU = {Qqu, Qbu} with Q^ as the 
L^ projection onto appropriately defined spaces. The following norms are used to 
measure the error in all of the numerical experiments: 

iJ^ semi-norm: |||eh||| = I /^ / jV^.e/ipdr + /i^^ / \Qbeo — eb\^ds 
Krp^j- Jt Jot 

Element-based L^ norm: ||eo|| = ( V" / \eo\'^dT 

7.1. On Triangular Mesh. Consider the second order elliptic equation that 
seeks an unknown function u = u(x, y) satisfying 

-V • (aV) = / 

in the square domain 51 = (0, 1) x (0, 1) with Dirichlet boundary condition. The 
boundary condition u|ao — g and / are chosen such that the exact solution is given 
by u = sin(7ra;) cos(7ry) and 

■ + y^ + 1 xy 

xy x^ + y"^ + 1 

The triangular mesh 7^ used in this example is constructed by: 1) uniformly 
partitioning the domain into nxn sub-rectangles; 2) dividing each rectangular element 
by the diagonal line with a negative slope. The mesh size is denoted by h = 1/n. The 
lowest order (fc = 1) weak Galerkin element is used for obtaining the weak Galerkin 
solution Uh — {uo,Ub}', i.e., uq and Ub are polynomials of degree k — 1 and degree 
fc — 1 = respectively on each element T G Th- 



Table 7.1 shows the convergence rate for WG solutions measured in H^ and L^ 
norms. The numerical results indicate that the WG solution of linear element is 
convergent with rate 0{h) in H^ and 0{h^) in L^ norms. 

In the second example, we consider the Poisson problem that seeks an unknown 
function u = u{x,y) satisfying 

-Au = f 

in the square domain fl — (0,1)^. Like the first example, the exact solution here is 
given by u = sin(7ra;) cos{TTy) and g and / are chosen accordingly to match the exact 
solution. 
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Table 7.1 
Example 1. Convergence rate of lowest order WG (k ■ 



1) on triangular meshes. 



h 


III eft III 


order 


Ik'oll 


order 


1/4 


1.3240e+00 




1.5784C+00 




1/8 


6.6333e-01 


9.9710e-01 


3.6890e-01 


2.0972 


1/16 


3.3182e-01 


9.9933e-01 


9.0622e-02 


2.0253 


1/32 


1.6593e-01 


9.9983e-01 


2.2556e-02 


2.0064 


1/64 


8.2966e-02 


9.9998e-01 


5.6326e-03 


2.0016 


1/128 


4.1483e-02 


1.0000 


1.4078e-03 


2.0004 



Table 7.2 
Example 2. Convergence rate of lowest order WG (k ■ 



1) on triangular meshes. 



h 


left III 


lleftll 


hhWsn 


1/2 


2.7935e-01 


6.1268e-01 


5.7099e-02 


1/4 


1.4354e-01 


1.5876e-01 


1.3892e-02 


1/8 


7.2436e-02 


4.0043C-02 


3.5430e-03 


1/16 


3.6315C-02 


1.0033e-02 


8.9325e-04 


1/32 


1.8170e-02 


2.5095e-03 


2.2384e-04 


1/64 


9.0865e-03 


6.2747e-04 


5.5994e-05 


1/128 


4.5435e-03 


1.5687C-04 


1.4001C-05 


0{h-),r^ 


9.9232C-01 


1.9913 


1.9955 



The very same triangular mesh is employed in the numerical calculation. As- 
sociated with this triangular mesh 7ft, two weak Galerkin elements with fc = 1 and 
k = 2 are used in the computation of the weak Galerkin finite element solution ut ■ For 
simplicity, these two elements shall be referred to as (Pi(T), Po(e)) and {P2{T) , Pi{e)) . 

Tables |7.2| and |7.3| show the numerical results on rate of convergence for the WG 
solutions in H^ and L? norms associated with k — 1 and k = 2, respectively. Note 
that ||e/j||£^ is a discrete L^ norm for the approximation Ub on the boundary of each 
element. Optimal rates of convergence are observed numerically for each case. 

7.2. On Quadrilateral Meshes. In this test, we solve the same poisson equa- 
tion considered in the second example by using quadrilateral meshes. We start with 



an initial quadrilateral mesh, shown as in Figure 7.1 (Left). The mesh is then suc- 
cessively refined by connecting the barycenter of each coarse element with the middle 
points of its edges, shown as in Figure 7.1 (Right). For the quadrilateral mesh 7ft, 



two weak Galerkin elements with fc = 1 and fc = 2 are used in the WG finite element 



scheme (3.4) 



Tables 7.4 and 7.5 show the rate of convergence for the WG solutions in H^ and 
L^ norms associated with fc = 1 and fc = 2 on quadrilateral meshes, respectively. 
Optimal rates of convergence are observed numerically. 

7.3. On Honeycomb Mesh. In the forth test, we solve the Poisson equation 
on the domain of unit square with exact solution u — sin(7ra::) sin(7ry). The Dirichlet 
boundary data g and / are chosen to match the exact solution. The numerical exper- 
iment is performed on the honeycomb mesh as shown in Figure |7.2[ The linear WG 
element (fc = 1) is used in this numerical computation. 



The error profile is presented in Table 7.6 which confirms the convergence rates 
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Table 7.3 
Example 2. Convergence rate of second order WG (k - 



■ 2) on triangular meshes. 



h 


lllehil 


W^hW 


hhWe,, 


1/2 


1.7886e-01 


9.4815e-02 


3.3742e-02 


1/4 


4.8010C-02 


1.2186e-02 


4.9969e-03 


1/8 


1.2327e-02 


1.5271e-03 


6.6539e-04 


1/16 


3.1139e-03 


1.9077e-04 


8.5226e-05 


1/32 


7.8188e-04 


2.3829e-05 


1.0763e-05 


1/64 


1.9586e-04 


2.9774C-06 


1.3516e-06 


1/128 


4.9009e-05 


3.7210C-07 


1.6932C-07 


0{h''),r = 


1.9769 


2.9956 


2.9453 



1 

0.8 








0.6 


~_\ 




0.4 


7 


\ 


\ 




0.2 






/ 




) 0.2 


0.4 


0.6 


0.8 1 




Fig. 7.1. Mesh level 1 (Left) and mesh level 2 (Right) for example 2. 

Table 7.4 
Example 3. Error and rate of convergence for first order WG on quadrilateral meshes. 



h 


\\K\\\ 


order 


lleoll 


order 


2.9350e-01 


1.9612C+00 




2.1072C+00 




1.4675e-01 


1.0349e+00 


9.2225e-01 


5.7219C-01 


1.8808 


7.3376e-02 


5.2434e-01 


9.8094e-01 


1.4458e-01 


1.9847 


3.6688e-02 


2.6323e-01 


9.9418e-01 


3.5655e-02 


2.0197 


1.8344e-02 


1.3179e-01 


9.9808C-01 


8.6047e-03 


2.0509 


9.1720e-03 


6.5925e-02 


9.9934e-01 


2.0184e-03 


2.0919 



predicted by the theory. 

7.4. On Deformed Cubic Meshes. In the fifth test, the Poisson equation is 
solved on a three dimensional domain ft — (0, 1)'^. The exact solution is chosen as 

u = sin(27ra:;) sin(27rj/) sin(27r2;), 

and the Dirichlet boundary date g and / are chosen accordingly to match the exact 
solution. 



Deformed cubic meshes are used in this test, see Figure 7.3 (Left) for an illustrative 



element. The construction of the deformed cubic mesh starts with a coarse mesh. The 
next level of mesh is derived by refining each deformed cube element into 8 sub-cubes, 
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Table 7.5 
Example 3. Error and rate of convergence for second order WG on quadrilateral meshes. 



h 


ll|e/i||| 


order 


l|eo|| 


order 


1/2 


1.7955e-01 




1.4891e-01 




1/4 


8.7059C-02 


1.0444 


1.8597C-02 


3.0013 


1/8 


2.8202C-02 


1.6262 


2.1311e-03 


3.1254 


1/16 


7.8114C-03 


1.8521 


2.4865e-04 


3.0995 


1/32 


2.0347e-03 


1.9408 


2.9964e-05 


3.0528 


1/64 


5.1767e-04 


1.9747 


3.6806e-06 


3.0252 


1/128 


1.3045C-04 


1.9885 


4.5627e-07 


3.0120 




Fig. 7.2. Honeycomb mesh for example 3. 

Table 7.6 
Example 4- Error and rate of convergence for linear WG element on honeycomb meshes. 



h 


lllehil 


order 


lleoll 


order 


1.6667e-01 


3.3201e-01 




1.6006e-02 




8.3333e-02 


1.6824e-01 


9.8067e-01 


3.9061e-03 


2.0347 


4.1667e-02 


8.4784e-02 


9.8867e-01 


9.6442e-04 


2.0180 


2.0833e-02 


4.2570e-02 


9.9392e-01 


2.3960e-04 


2.0090 


1.0417e-02 


2.1331e-02 


9.9695e-01 


5.9711e-05 


2.0047 


5.2083e-03 


1.0677e-02 


9.9839e-01 


1.4904e-05 


2.0022 



as shown in Figure |7.3| (Right) . Table 7.7 reports some numerical results for different 
level of meshes. It can be seen that a convergent rate of 0{h) in H^ and 0(h^) in L"^ 
norms are achieved for the corresponding WG finite element solutions. This confirms 
the theory developed in earlier sections. 
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-1 -1 
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FlG. 7.3. Mesh level 1 (Left) and mesh level 2 (Right) for example 4- 

Table 7.7 
Example 5. Error and convergence rate for k = 1 on deformed cubic mesh. 



h 


lllehil 


order 


l|eo|| 


order 


1/2 


5.7522 




9.1990 




1/4 


1.3332 


2.1092 


1.5684 


2.5522 


1/8 


6.4071e-01 


1.0571 


2.7495e-01 


2.5121 


1/16 


3.2398e-01 


9.8377e-01 


6.8687e-02 


2.0011 


1/32 


1.6201e-01 


9.9982C-01 


1.7150e-02 


2.0018 
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